* Replication of Figure 6, panel a, using Golder marginal effects plot
* You will need to specify the directory in which the repliaction data are stored in the second line of code 

clear
cd "[Directory with replication data here]"

version 11.0
#delimit ;
log using "sim6a.log", replace;
set more off;

use "fordham_wp2025.dta", clear ;

logit multilat black time interactb ag interacta lynchings interactl cotton textapp, cluster(icpsr);

preserve;

set seed 339487731;

drawnorm SN_b1-SN_b10, n(10000) means(e(b)) cov(e(V)) clear;

postutil clear;
postfile mypost prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi 
            using "sim6a.dta", replace;
            
noisily display "start";

local a=0 ;
while `a' <= 36 { ;

{;

scalar h_time=17;
scalar h_cotton=0.26;
scalar h_textapp=4.7;
scalar h_black=20;
scalar h_ag=20;
scalar h_lynchings=30;
scalar h_constant=1;

    generate x_betahat0  = SN_b1*h_black 
                            + SN_b2*h_time  
                            + SN_b3*h_black*h_time
                            + SN_b4*h_ag 
                            + SN_b5*h_ag*h_time
                            + SN_b6*h_lynchings 
                            + SN_b7*h_lynchings*h_time
                            + SN_b8*h_cotton
                            + SN_b9*`a'
                            + SN_b10*h_constant;
                            
    generate x_betahat1  = SN_b1*h_black
                            + SN_b2*(h_time+6)
                            + SN_b3*h_black*(h_time+6) 
                            + SN_b4*h_ag 
                            + SN_b5*h_ag*(h_time+6) 
                            + SN_b6*h_lynchings 
                            + SN_b7*h_lynchings*(h_time+6) 
                            + SN_b8*h_cotton
                            + SN_b9*`a' 
                            + SN_b10*h_constant;
    
    gen prob0 =normal(x_betahat0);
    gen prob1=normal(x_betahat1);
    gen diff=prob1-prob0;
    
    egen probhat0 =mean(prob0);
    egen probhat1=mean(prob1);
    egen diffhat=mean(diff);
    
    tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi;  

    _pctile prob0, p(2.5,97.5) ;
    scalar `lo0' = r(r1);
    scalar `hi0' = r(r2);  
    
    _pctile prob1, p(2.5,97.5);
    scalar `lo1'= r(r1);
    scalar `hi1'= r(r2);  
    
    _pctile diff, p(2.5,97.5);
    scalar `diff_lo'= r(r1);
    scalar `diff_hi'= r(r2);  
   
    scalar `prob_hat0'=probhat0;
    scalar `prob_hat1'=probhat1;
    scalar `diff_hat'=diffhat;
    
    post mypost (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') (`lo1') (`hi1') 
                (`diff_hat') (`diff_lo') (`diff_hi');
   
    };      
        
    drop x_betahat0 x_betahat1 prob0 prob1 diff probhat0 probhat1 diffhat;
	
    local a=`a'+ 1;
    
    display "." _c;
    
} ;

display "";

postclose mypost;

restore;

merge using "sim6a.dta";

gen MV = _n-1;

replace  MV=. if _n>37;

graph twoway hist textapp, percent color(gs14) yaxis(2)
        ||  line prob_hat0 MV, clwidth(medium) clcolor(blue) clcolor(black)
        ||  line lo0 MV, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line hi0 MV, clpattern(dash) clwidth(thin) clcolor(black)
        ||  ,   
            xlabel(0 10 20 30 36, nogrid labsize(3)) 
            ylabel(0 0.2 0.4 0.6 0.8 1, axis(1) nogrid labsize(3))
            ylabel(0 5 10 15 20 25 30 35, axis(2) nogrid labsize(3))
            yscale(noline alt) yscale(noline alt axis(2)) xscale(noline) legend(off)
            yline(0, lcolor(black))
            yline(0 0.2 0.4 0.6 0.8 1, lcolor(white))
            xtitle("Percent of Workforce in Textiles and Apparel Sectors", size(3))
            ytitle("Probability of Voting in Support of Multilateralism", size(3))
            ytitle("Percent of Observations" , axis(2) size(3))
            xsca(titlegap(4)) ysca(titlegap(4)) ysca(axis(2) titlegap(4))
            scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));
            
log close;
